% Copyright (c) 2012, John Gebbie. All rights reserved.
%
% This library is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation; either version 2.1 of the License, or (at
% your option) any later version.
%
% This library is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
% General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this library; if not, write to the Free Software Foundation,
% Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA


classdef ortrace < handle
    
    properties (SetAccess = private)
        
        %------------------------------------------------------------------
        % Environment parameters
        %------------------------------------------------------------------
        
        params
        
        %------------------------------------------------------------------
        % SSP properties
        %------------------------------------------------------------------
        
        ssp_c
        ssp_z
        ssp_r
        ssp_dcdz
        ssp_step_z
        
        %------------------------------------------------------------------
        % Output
        %------------------------------------------------------------------
        
        rays_out
        ray_out
        
    end % properties
    
    
    properties (Constant)
        
        POS_R = 1;      % range
        POS_Z = 2;      % depth
        POS_XI = 3;     % xi - drds/c
        POS_ZETA = 4;   % zeta - dzds/c
        POS_DTRV = 5;   % distance traveled (meters)
        POS_TTRV = 6;   % time traveled (seconds)
        POS_RC = 7;     % accumulated reflection coefficient
        POS_NBS = 8;    % number of surface bounces
        POS_NBB = 9;    % number of bottom bounces
        POS_RIDX = 10;  % index in output range vector
        POS_BNCS = 11;  % entry is a surface bounce
        POS_BNCB = 12;  % entry is a bottom bounce
        POS_GANG = 13;  % grazing angle for bounce (NaN if not a bounce)
        POS_SIDX = 14;  % seabed index (optimization)
        NPOS = 14;
        
    end % properties
    
    
    methods
        
        function init_environment(obj, params)
            
            obj.params = params;
            
            obj.discretize_ssp();
            
        end % function
        
        
        function [] = cast_ray_fan(obj)
            
            if isfield(obj.params, 'thetas')
                thetas = obj.params.thetas;
                theta_n = length(thetas);
            else
                theta_mn = obj.params.theta_mn_mx_n(1);
                theta_mx = obj.params.theta_mn_mx_n(2);
                theta_n = obj.params.theta_mn_mx_n(3);
                thetas = linspace(theta_mn,theta_mx,theta_n);
            end
            
            obj.rays_out = cell(theta_n,1);
            
            for ix = 1:length(thetas)
                obj.cast_ray(thetas(ix));
                obj.rays_out{ix} = obj.ray_out;
            end
            
        end % function
        
        
        function [] = cast_ray(obj, theta_0)
            
            %--------------------------------------------------------------
            % Input parameters
            %--------------------------------------------------------------
            
            rs = obj.params.rs;
            zs = obj.params.zs;
            orng_mn = obj.params.output_ranges_mn_mx_n(1);
            orng_mx = obj.params.output_ranges_mn_mx_n(2);
            orng_n = obj.params.output_ranges_mn_mx_n(3);
            if orng_n > 1
                orng_seg = (orng_mx-orng_mn)/(orng_n-1);
            else
                orng_seg = 1;
            end
            output_full_ray = obj.params.output_full_ray;
            output_grazing_angle = obj.params.output_grazing_angle;
            step_m = obj.params.step_m;
            min_amplitude_threshold = 10^(-obj.params.max_db_loss/20);
            max_depth = obj.params.max_depth;
            
            %--------------------------------------------------------------
            % Init ray data struct
            %--------------------------------------------------------------
            
            ray = nan(ortrace.NPOS, 1);
            ray(ortrace.POS_R)    = rs;
            ray(ortrace.POS_Z)    = zs;
            tmp_cs = interp1(obj.ssp_z, obj.ssp_c(:, 1), zs);
            ray(ortrace.POS_XI)   = cos(theta_0)./tmp_cs;
            ray(ortrace.POS_ZETA) = sin(theta_0)./tmp_cs;
            ray(ortrace.POS_DTRV) = 0;
            ray(ortrace.POS_TTRV) = 0;
            ray(ortrace.POS_RC)   = 1;
            ray(ortrace.POS_NBS)  = 0;
            ray(ortrace.POS_NBB)  = 0;
            ray(ortrace.POS_RIDX) = NaN;
            ray(ortrace.POS_BNCS) = NaN;
            ray(ortrace.POS_BNCB) = NaN;
            ray(ortrace.POS_GANG) = NaN;
            ray(ortrace.POS_SIDX) = 1;
            
            obj.ray_out = [];
            
            if output_full_ray
                obj.ray_out(:, end+1) = ray;
            end
            
            %--------------------------------------------------------------
            % Start propagating ray forward
            %--------------------------------------------------------------
            
            range_idx_b = 1;
            range_idx_a = 1;
            
            while ...
                    ... % ray is in valid positive range
                    (0 <= ray(ortrace.POS_R) && ...
                    ray(ortrace.POS_R) <= orng_mx)...
                    && ... % ray has not attenuated past threshold
                    ((ray(ortrace.POS_DTRV) <= 0) || ...
                    (abs(ray(ortrace.POS_RC)) / ...
                    ray(ortrace.POS_DTRV) >= ...
                    min_amplitude_threshold)) ...
                    && ... % ray is in valid depth region
                    (ray(ortrace.POS_Z) <= max_depth)
                
                %----------------------------------------------------------
                % Current depth converted to an integer to be used
                %   for indexing
                %----------------------------------------------------------
                depth_idx = ...
                    floor(ray(ortrace.POS_Z)/obj.ssp_step_z) + 1;
                depth_idx = max(depth_idx, 1);
                depth_idx = min(depth_idx, length(obj.ssp_c));
                
                %----------------------------------------------------------
                % Find ranges in ssp that ray is currently between
                %----------------------------------------------------------
                while true
                    if ray(ortrace.POS_R) < obj.ssp_r(range_idx_b)
                        break;
                    end
                    if range_idx_b == length(obj.ssp_r)
                        range_idx_a = range_idx_b;
                        break;
                    end
                    range_idx_b = range_idx_b + 1;
                    range_idx_a = range_idx_b - 1;
                end
                
                %----------------------------------------------------------
                % Sound speed and derivative at current depth
                %----------------------------------------------------------
                if range_idx_a == range_idx_b
                    c = obj.ssp_c(depth_idx, range_idx_a);
                    dcdz = obj.ssp_dcdz(depth_idx, range_idx_a);
                    dcdr = 0;
                else
                    tmp_ra = obj.ssp_r(range_idx_a);
                    tmp_rb = obj.ssp_r(range_idx_b);
                    % fraction of distance between a and b
                    tmp_t = (ray(ortrace.POS_R)-tmp_ra)/(tmp_rb-tmp_ra);
                    c = obj.ssp_c(depth_idx, [range_idx_a range_idx_b]);
                    dcdr = (c(2)-c(1))/(tmp_rb-tmp_ra);
                    c = c(1).*(1-tmp_t) + c(2).*tmp_t;
                    dcdz = obj.ssp_dcdz(depth_idx, ...
                        [range_idx_a range_idx_b]);
                    dcdz = dcdz(1).*(1-tmp_t) + dcdz(2).*tmp_t;
                end
                
                %----------------------------------------------------------
                % Derivative with respect to 's'
                %   - 's' is distance along the ray
                %----------------------------------------------------------
                x = nan(ortrace.NPOS, 1);
                x(ortrace.POS_R)    = c*ray(ortrace.POS_XI);
                x(ortrace.POS_Z)    = c*ray(ortrace.POS_ZETA);
                x(ortrace.POS_XI)   = -dcdr/(c.^2);
                x(ortrace.POS_ZETA) = -dcdz/(c.^2);
                x(ortrace.POS_DTRV) = 1;
                x(ortrace.POS_TTRV) = 1./c;
                x(ortrace.POS_RC)   = 0;
                x(ortrace.POS_NBS)  = 0;
                x(ortrace.POS_NBB)  = 0;
                x(ortrace.POS_RIDX) = NaN;
                x(ortrace.POS_BNCS) = NaN;
                x(ortrace.POS_BNCB) = NaN;
                x(ortrace.POS_GANG) = NaN;
                x(ortrace.POS_SIDX) = 0;
                
                %----------------------------------------------------------
                % Step size (in meters)
                %
                % If ray segment crosses or ends on an output grid range
                % boundary, only propagate the ray to this point and outut
                % the ray state.
                %
                % The range of the k'th output grid point is
                %   (k-1)*(mx-mn)/(n-1)+mn              (n > 1) mn
                %   (n == 1)
                %
                % The strategy is to compute the output grid point
                % immediately to the left of the start and end points of
                % the ray. If they differ, output to the range of the
                % output grid point.
                %
                % At range 'r', the index of the output grid point
                % immediately to the left is given as
                %   floor(r/mn)                         (n == 1)
                %   floor((r-mn)/((mx-mn)/(n-1)))+1     (n > 1)
                % If the resulting index is between 1 and n, the ray ends
                % on an output grid point and is output, otherwise ignored.
                %
                % To handle outputting to a grid point, the range of the
                % grid point is
                %   mn                                  (n == 1)
                %   mn+(i-1)*(mx-mn)/(n-1)              (n > 1)
                % where 'i' is the index of the start or end of the ray,
                % depending on which ever falls between two. The new 'h'
                % (ray segment length) is
                %   abs((g-r1)/x)
                %----------------------------------------------------------
                
                h = step_m;
                at_output_range = false;
                output_range_idx = NaN;
                if orng_n > 0
                    tmp_r1 = ray(ortrace.POS_R);
                    tmp_r2 = tmp_r1 + h * x(ortrace.POS_R);
                    
                    tmp_gate = @(X,MN,MX) max(min(X,MX),MN);
                    % add or subtract 1e-9 to avoid rounding errors which
                    % may cause the start of the ray on next iteration to
                    % map to the previous index.
                    if x(ortrace.POS_R) >= 0
                        tmp_opn = @(X) tmp_gate(floor(X+1e-9)+1,0,orng_n);
                        tmp_opi = @(X) X+1;
                    else
                        tmp_opn = @(X) tmp_gate(ceil(X-1e-9),0,orng_n);
                        tmp_opi = @(X) X;
                    end
                    
                    tmp_i1 = tmp_opn((tmp_r1-orng_mn) / orng_seg);
                    tmp_i2 = tmp_opn((tmp_r2-orng_mn) / orng_seg);
                    
                    at_output_range = tmp_i1 ~= tmp_i2;
                    
                    if at_output_range
                        if orng_n == 1
                            output_range_idx = 1;
                            tmp_g = orng_mn;
                        else
                            output_range_idx = tmp_opi(tmp_i1);
                            tmp_g = orng_mn + ...
                                (output_range_idx-1)*...
                                (orng_mx-orng_mn)/(orng_n-1);
                        end
                        h = abs( (tmp_g - tmp_r1) / x(ortrace.POS_R) );
                        assert(h>0);
                    end
                end
                
                %----------------------------------------------------------
                % Step forward
                %----------------------------------------------------------
                ray_next = ray + h*x;
                
                %----------------------------------------------------------
                % Check for seabed or surface boundary intersection
                %----------------------------------------------------------
                [b_st, ray_next] = get_boundary_intersection(obj, ...
                    ray, ray_next, x);
                
                %----------------------------------------------------------
                % Boundary reflection
                %----------------------------------------------------------
                if b_st.at_boundary
                    
                    % ray_next is point at which the ray interects the
                    % boundary (it is colinear with boundary segment),
                    % which we call point 'b'. To find the reflected angle,
                    % we first define the boundary segment starting at R/Z
                    % of point 'b' and extending to point 'c' parallel to
                    % the boundary (length does not matter). The angle of
                    % the incident ray is determined by xi and zeta at
                    % ray_next, so we compute the start of the incident ray
                    % by subtracting [xi;zeta] from point 'b' to point 'a'.
                    % We are now set up for using our ray reflection
                    % routine.
                    
                    %------------------------------------------------------
                    % Reflect ray (compute new direction which is defined
                    % by xi and zeta)
                    %------------------------------------------------------
                    ax = -ray_next(ortrace.POS_XI);
                    ay = -ray_next(ortrace.POS_ZETA);
                    bx = b_st.boundary_slope_dr;
                    by = b_st.boundary_slope_dz;
                    refl_xy = ortrace.ray_reflect_from_line(...
                        ax, ay, bx, by);
                    ray_next(ortrace.POS_XI)   = refl_xy(1);
                    ray_next(ortrace.POS_ZETA) = refl_xy(2);
                    
                    %------------------------------------------------------
                    % Compute angle of incidence (measured with respect
                    %   to boundary)
                    %------------------------------------------------------
                    % take dot product of incident and reflected rays which
                    % gives the cosine of the angle between them. since
                    % angle of incidence equals the angle of reflection,
                    % divide by 2 to get angle of incidence with respect to
                    % the normal. subtract from pi/2 to get angle with
                    % respect to boundary.
                    v1 = [ax ay];
                    v1 = v1 ./ sqrt(sum(v1.^2));
                    v2 = refl_xy;
                    v2 = v2 ./ sqrt(sum(v2.^2));
                    graze_ang = pi/2-acos(v1(1)*v2(1)+v1(2)*v2(2))/2;
                    
                    %------------------------------------------------------
                    % Check for total internal reflection (skip if incident
                    % angle is shallow enough since all energy is
                    % reflected)
                    %------------------------------------------------------
                    cos_theta1 = b_st.boundary_c/c * cos(graze_ang);
                    if cos_theta1 <= 1
                        
                        %--------------------------------------------------
                        % Compute angle of transmission (measured with
                        %   respect to horizontal)
                        %--------------------------------------------------
                        theta2 = acos(b_st.boundary_c/c * cos(graze_ang));
                        
                        %--------------------------------------------------
                        % Compute the impedance of water and seabed
                        %   for the given angles
                        %--------------------------------------------------
                        Z2 = b_st.boundary_impedance / sin(theta2);
                        Z1 = 1.0 * c / sin(graze_ang);
                        
                        %--------------------------------------------------
                        % Accmululate the coefficient in the 6th
                        %   vector slot
                        %--------------------------------------------------
                        R = (Z2 - Z1) / (Z2 + Z1);
                        ray_next(ortrace.POS_RC) = ...
                            ray(ortrace.POS_RC) * R;
                    end
                    
                    %------------------------------------------------------
                    % If previously propagating ray to output range grid
                    % point, skip this now that boundary was intersected.
                    %------------------------------------------------------
                    at_output_range = false;
                    output_range_idx = NaN;
                    
                else
                    
                    %------------------------------------------------------
                    % No boundary intersection, so no grazing angle to
                    % output
                    %------------------------------------------------------
                    graze_ang = NaN;
                end
                
                %----------------------------------------------------------
                % Write some additional output variables
                %----------------------------------------------------------
                ray_next(ortrace.POS_RIDX) = output_range_idx;
                ray_next(ortrace.POS_GANG) = graze_ang;
                ray_next(ortrace.POS_BNCS) = ...
                    ray(ortrace.POS_NBS) ~= ray_next(ortrace.POS_NBS);
                ray_next(ortrace.POS_BNCB) = ...
                    ray(ortrace.POS_NBB) ~= ray_next(ortrace.POS_NBB);
                
                %----------------------------------------------------------
                % Store the next step in the ray at the end of the array
                %----------------------------------------------------------
                ray = ray_next;
                if at_output_range || ...
                        output_full_ray || ...
                        (output_grazing_angle && ~isnan(graze_ang))
                    obj.ray_out(:, end+1) = ray_next;
                end
            end
            
        end % function
        
        
        function [arrivals] = get_arrivals(obj, r_idx, z)
            
            % Define a "type" of ray has having a given number of surface
            % and bottom bounces.
            %
            % Adjacent rays of the same type can be combined by linear
            % interpolation to estimate the arrival information at any
            % point falling between those rays. Interpolated quantities are
            % distance traveled and time traveled.  Grazing angles can also
            % be interpolated, but since multiple boundary intersections
            % can exist for each pair of adjacent rays, each must be
            % interpolated separately.
            %
            % If the requested depth falls between two adjacent rays of
            % different types,
            %
            % If the requested depth falls beyond the first or last ray,
            % output nothing, or an error.
            
        end % function
        
    end % methods
    
    
    methods (Access = private)
        
        function [b_struct, ray_next] = get_boundary_intersection(obj, ...
                ray, ray_next, x)
            
            %--------------------------------------------------------------
            % Check for bottom intersection
            %
            %   Test each segment in the seabed vector having
            %   overlapping ranges with the range span of the ray.
            %--------------------------------------------------------------
            
            %--------------------------------------------------------------
            % Get the indices in the seabed vector that the ray
            % falls between.
            %--------------------------------------------------------------
            [idx_seabed_a, idx_seabed_b] = ...
                obj.find_seabed_indices(...
                ray(ortrace.POS_R), ...
                ray_next(ortrace.POS_R), ...
                ray(ortrace.POS_SIDX));
            
            %--------------------------------------------------------------
            % Save the last seabed index for next iteration
            %--------------------------------------------------------------
            ray(ortrace.POS_SIDX) = max(1, idx_seabed_a);
            
            %--------------------------------------------------------------
            % Setup line-segment for ray path
            %--------------------------------------------------------------
            seg_cr = ray(ortrace.POS_R);
            seg_cz = ray(ortrace.POS_Z);
            seg_dr = ray_next(ortrace.POS_R);
            seg_dz = ray_next(ortrace.POS_Z);
            
            %--------------------------------------------------------------
            % Setup the direction over which to check seabed
            % segments.  For positively traveling rays, search in
            % the positive direction, and vice-versa.
            %--------------------------------------------------------------
            if x(ortrace.POS_R) >= 0
                idx_seabed_start = idx_seabed_a;
                idx_seabed_end = idx_seabed_b-1;
                idx_direction = 1;
            else
                idx_seabed_start = idx_seabed_b-1;
                idx_seabed_end = idx_seabed_a;
                idx_direction = -1;
            end
            
            %--------------------------------------------------------------
            % Loop through all seabed segments looking for
            % intersection with ray
            %--------------------------------------------------------------
            found_intersection = false;
            for idx_seabed = idx_seabed_start:idx_direction:idx_seabed_end
                
                %----------------------------------------------------------
                % Create special seabed segments when the ray lies
                % past the end of defined region of seabed
                %----------------------------------------------------------
                if length(obj.params.seabed_range) <= idx_seabed
                    boundary_seg_ar = obj.params.seabed_range(end);
                    boundary_seg_az = obj.params.seabed_depth(end);
                    boundary_seg_br = max([seg_cr, seg_dr])+1;
                    boundary_seg_bz = boundary_seg_az;
                    idx_seabed_a = length(obj.params.seabed_range);
                    idx_seabed_b = length(obj.params.seabed_range);
                elseif idx_seabed < 1
                    boundary_seg_ar = obj.params.seabed_range(1);
                    boundary_seg_az = obj.params.seabed_depth(1);
                    boundary_seg_br = min([seg_cr, seg_dr])-1;
                    boundary_seg_bz = boundary_seg_az;
                    idx_seabed_a = 1;
                    idx_seabed_b = 1;
                else
                    boundary_seg_ar = obj.params.seabed_range(idx_seabed);
                    boundary_seg_az = obj.params.seabed_depth(idx_seabed);
                    boundary_seg_br = obj.params.seabed_range(idx_seabed+1);
                    boundary_seg_bz = obj.params.seabed_depth(idx_seabed+1);
                    idx_seabed_a = idx_seabed;
                    idx_seabed_b = idx_seabed+1;
                end
                
                if ~ortrace.colinear(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_cr, seg_cz) ...
                        && ...
                        ortrace.line_segments_intersect(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_cr, seg_cz, seg_dr, seg_dz) ...
                        || ...
                        ortrace.colinear(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_dr, seg_dz) ...
                        && ... % extend segment and check again for i/sect
                        ortrace.line_segments_intersect(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_cr, seg_cz, 2*seg_dr-seg_cr, 2*seg_dz-seg_cz)
                    
                    found_intersection = true;
                    break;
                end
                
            end
            
            %--------------------------------------------------------------
            % If intersection with seabed is found, propagate ray
            % to seabed
            %--------------------------------------------------------------
            if found_intersection
                
                %----------------------------------------------------------
                % Get intersection point
                %----------------------------------------------------------
                rz_end = ...
                    ortrace.line_line_intersection_point(...
                    boundary_seg_ar, boundary_seg_az, ...
                    boundary_seg_br, boundary_seg_bz, ...
                    seg_cr, seg_cz, seg_dr, seg_dz);
                rz_start = ray([ortrace.POS_R, ortrace.POS_Z]);
                
                %----------------------------------------------------------
                % Compute new step size from ray to tmp_xz
                %----------------------------------------------------------
                dist = sqrt(sum((rz_end-rz_start).^2));
                
                %----------------------------------------------------------
                % Step forward to boundary
                %----------------------------------------------------------
                ray_next = ray + dist*x;
                % avoid rounding errors..
                ray_next(ortrace.POS_R) = rz_end(1);
                ray_next(ortrace.POS_Z) = rz_end(2);
                
                %----------------------------------------------------------
                % Record bottom intersection
                %----------------------------------------------------------
                ray_next(ortrace.POS_NBB) = ray(ortrace.POS_NBB) + 1;
                
                %----------------------------------------------------------
                % Setup for boundary reflection (common code below)
                %----------------------------------------------------------
                at_boundary = true;
                if idx_seabed_a == idx_seabed_b
                    boundary_c = obj.params.seabed_c(idx_seabed_a);
                    boundary_rho = obj.params.seabed_rho(idx_seabed_a);
                    boundary_impedance = boundary_rho*boundary_c;
                    boundary_slope_dr = 1;
                    boundary_slope_dz = 0;
                else
                    tmp_a = obj.params.seabed_range(idx_seabed_b) - ...
                        obj.params.seabed_range(idx_seabed_a);
                    tmp_b = ray_next(ortrace.POS_R) - ...
                        obj.params.seabed_range(idx_seabed_a);
                    tmp_t = tmp_b / tmp_a;
                    boundary_c = ...
                        (1-tmp_t)*obj.params.seabed_c(idx_seabed_a) + ...
                        tmp_t*    obj.params.seabed_c(idx_seabed_b);
                    boundary_rho = ...
                        (1-tmp_t)*obj.params.seabed_rho(idx_seabed_a) + ...
                        tmp_t*    obj.params.seabed_rho(idx_seabed_b);
                    boundary_impedance = boundary_rho*boundary_c;
                    boundary_slope_dr = ...
                        obj.params.seabed_range(idx_seabed_b) - ...
                        obj.params.seabed_range(idx_seabed_a);
                    boundary_slope_dz = ...
                        obj.params.seabed_depth(idx_seabed_b) - ...
                        obj.params.seabed_depth(idx_seabed_a);
                end
                
            else
                
                %----------------------------------------------------------
                % Check for surface intersection
                %----------------------------------------------------------
                at_boundary = false;
                boundary_seg_ar = 0;
                boundary_seg_az = 0;
                boundary_seg_br = 1;
                boundary_seg_bz = 0;
                seg_cr = ray(ortrace.POS_R);
                seg_cz = ray(ortrace.POS_Z);
                seg_dr = ray_next(ortrace.POS_R);
                seg_dz = ray_next(ortrace.POS_Z);
                if ~ortrace.colinear(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_cr, seg_cz) ...
                        && ...
                        ~ortrace.line_segment_same_halfspace(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_cr, seg_cz, ...
                        seg_dr, seg_dz) ...
                        || ...
                        ortrace.colinear(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_dr, seg_dz)
                    
                    %------------------------------------------------------
                    % Get intersection point
                    %------------------------------------------------------
                    rz_end = ortrace.line_line_intersection_point(...
                        boundary_seg_ar, boundary_seg_az, ...
                        boundary_seg_br, boundary_seg_bz, ...
                        seg_cr, seg_cz, seg_dr, seg_dz);
                    rz_start = ray([ortrace.POS_R, ortrace.POS_Z]);
                    
                    %------------------------------------------------------
                    % Compute new step size from ray to tmp_xz
                    %------------------------------------------------------
                    dist = sqrt(sum((rz_end-rz_start).^2));
                    
                    %------------------------------------------------------
                    % Step forward to boundary
                    %------------------------------------------------------
                    ray_next = ray + dist*x;
                    % avoid rounding errors..
                    ray_next(ortrace.POS_R) = rz_end(1);
                    ray_next(ortrace.POS_Z) = rz_end(2);
                    
                    %------------------------------------------------------
                    % Record surface intersection
                    %------------------------------------------------------
                    ray_next(ortrace.POS_NBS) = ray(ortrace.POS_NBS) + 1;
                    
                    %------------------------------------------------------
                    % Setup for boundary reflection
                    %------------------------------------------------------
                    at_boundary = true;
                    boundary_c = 0;
                    boundary_impedance = 0; % pressure-release
                    boundary_slope_dr = 1;
                    boundary_slope_dz = 0;
                    
                end
            end
            
            b_struct.at_boundary                = at_boundary;
            if at_boundary
                b_struct.boundary_c             = boundary_c;
                b_struct.boundary_impedance     = boundary_impedance;
                b_struct.boundary_slope_dr      = boundary_slope_dr;
                b_struct.boundary_slope_dz      = boundary_slope_dz;
            end
            
        end % function
        
        
        function discretize_ssp(obj)
            
            %--------------------------------------------------------------
            % Sort the input SSP
            %--------------------------------------------------------------
            
            [issp_r, tmp_range_idx] = sort(obj.params.ssp_r);
            issp_cz = obj.params.ssp_cz(tmp_range_idx);
            for iRng = 1:length(issp_cz)
                [issp_cz{iRng}.z, tmp_depths_idx] = sort(issp_cz{iRng}.z);
                issp_cz{iRng}.c = issp_cz{iRng}.c(tmp_depths_idx);
            end
            
            %--------------------------------------------------------------
            % Discretize SSP
            %--------------------------------------------------------------
            
            % save the range vector
            obj.ssp_r = issp_r;
            % save the step size
            obj.ssp_step_z = obj.params.ssp_step_z;
            % depths every ssp_step_z meters
            obj.ssp_z = (0 : obj.ssp_step_z : ...
                ceil(max(obj.params.seabed_depth)/obj.ssp_step_z) * ...
                obj.ssp_step_z).';
            obj.ssp_c  = nan(length(obj.ssp_z), length(issp_r));
            % if the deepest sound speed in the input ssp is less than the
            % deepest sound speed in the output buffer, artifically extend
            % the sound speed from the input sound speed to the deepest
            % level of the output. this ensures that each depth in the
            % input profile falls between two points in the input profile.
            for iRng = 1:length(issp_r)
                if issp_cz{iRng}.z(end) < obj.ssp_z(end)
                    issp_cz{iRng}.c(end+1) = issp_cz{iRng}.c(end);
                    issp_cz{iRng}.z(end+1) = obj.ssp_z(end);
                end
            end
            % for each output sound speed depth, search through the input
            % profile to find the adjacent pair of depths that the output
            % sound speed depth falls between, then linearly interpolate
            % between them. Assume the input sound speed profile is sorted.
            for iRng = 1:length(issp_r)
                iInp = 1;
                for iOut = 1:length(obj.ssp_z)
                    % search for pair of rows in the output ssp depths that
                    %   the current input ssp depth value falls between
                    while true
                        if issp_cz{iRng}.z(iInp) <= obj.ssp_z(iOut) && ...
                                obj.ssp_z(iOut) < issp_cz{iRng}.z(iInp+1)
                            break;
                        end
                        if length(issp_cz{iRng}.z)-1 <= iInp
                            break;
                        end
                        iInp = iInp + 1;
                    end
                    % linearly interpolate between the two depth points in
                    % the input ssp
                    tmp_z12 = issp_cz{iRng}.z(iInp:(iInp+1));
                    tmp_c12 = issp_cz{iRng}.c(iInp:(iInp+1));
                    tmp_slope = diff(tmp_c12)/diff(tmp_z12);
                    tmp_zoffset = obj.ssp_z(iOut) - tmp_z12(1);
                    obj.ssp_c(iOut, iRng) = ...
                        tmp_c12(1) + tmp_zoffset * tmp_slope;
                end
            end
            
            %--------------------------------------------------------------
            % Derivative of SSP - over depth
            %--------------------------------------------------------------
            
            obj.ssp_dcdz = diff(obj.ssp_c, [], 1) ./ obj.ssp_step_z;
            obj.ssp_dcdz = [...
                obj.ssp_dcdz(1, :); ...
                (obj.ssp_dcdz(1:end-1, :) + obj.ssp_dcdz(2:end, :))/2; ...
                obj.ssp_dcdz(end, :)];
            if size(obj.ssp_dcdz, 1) > 2
                % TODO: what does this do?
                obj.ssp_dcdz([end-1 end], :) = ...
                    repmat(obj.ssp_dcdz(end-2, :), 2, 1);
            end
            
        end % function
        
        
        function [ia, ib] = find_seabed_indices(obj, rng1, rng2, ia)
            % walk backward then forward in seabed vector to find the set
            % of indices in obj.seabed_range that are less than rng1 and
            % greater than rng2.  idx values are not advanced past the end
            % of the seabed_range vector.
            
            % swap, if needed
            if rng2 < rng1
                tmp = rng2;
                rng2 = rng1;
                rng1 = tmp;
            end
            
            n = length(obj.params.seabed_range);
            
            % walk 'a' backward until it is less than rng1
            while (1 < ia) && ~(obj.params.seabed_range(ia) <= rng1)
                ia = ia - 1;
            end
            
            % walk 'a' forward until another step would put it greater than
            % or equal to rng1
            while (ia < n) && ~(rng1 < obj.params.seabed_range(ia+1))
                ia = ia + 1;
            end
            
            ib = ia;
            
            % walk 'b' forward until it is greater than rng2
            while (ib < n) && ~(rng2 < obj.params.seabed_range(ib))
                ib = ib + 1;
            end
            
            % use dummy indices if range is outside of defined seabed
            % vector
            if rng1 < obj.params.seabed_range(ia)
                ia = ia - 1;
            end
            if obj.params.seabed_range(ib) < rng2
                ib = ib + 1;
            end
        end % function
        
    end % methods
    
    
    methods (Static, Access = private)
        
        function [retval] = ray_reflect_from_line(ax, ay, bx, by)
            % Returns ray reflection from a line. Reflection point is the
            % origin. Incident ray originates at point 'a' and intersects
            % line at the origin. The line is defined by point 'b' and the
            % origin. Point 'r' such that the origin to point 'r' is the
            % reflected ray and is of same length as the incident line
            % segment.
            
            % compute sin and cos of theta1 (relative to positive x axis)
            len = sqrt(bx^2+by^2);
            cos_theta1 = bx/len;
            sin_theta1 = by/len;
            % compute sin and cos of -theta1
            cos_mtheta1 = cos_theta1;
            sin_mtheta1 = -sin_theta1;
            % make rotation matrices for theta1 and -theta1 as well as
            % matrix to reflect over y-axis
            Mrotneg = [ cos_mtheta1 -sin_mtheta1 ; ...
                sin_mtheta1 cos_mtheta1 ];
            Mrefly  = [ -1 0 ; 0 1 ];
            Mrotpos = [ cos_theta1  -sin_theta1  ; ...
                sin_theta1  cos_theta1 ];
            % rotate bc by -theta1 flip over y axis rotate bc again by
            % +theta1
            retval = Mrotpos*Mrefly*Mrotneg*[ax ; ay];
        end % function
        
        
        function [retval] = line_line_intersection_point(...
                ax, ay, bx, by, cx, cy, dx, dy)
            % Returns point at which two lines intersect.  Lines are
            % defined by ab and cd.  NOTE: does not handle degenerate cases
            % of colinearity or colocation.
            %
            % Solve system of equations for each line.
            M = [   (by-ay),    -(bx-ax) ; ...
                (dy-cy),    -(dx-cx) ];
            b = [ ax*(by-ay)-ay*(bx-ax) ; cx*(dy-cy)-cy*(dx-cx) ];
            retval = M\b;
        end % function
        
        
        function [retval] = colinear(ax, ay, bx, by, cx, cy)
            % Returns true if three points are colinear
            %
            % cross-product |ab X ac| == 0
            abx = bx - ax;
            aby = by - ay;
            acx = cx - ax;
            acy = cy - ay;
            retval = abs(abx*acy-aby*acx) < ...
                (sqrt(abx^2+aby^2) + sqrt(acx^2+acy^2))/2 * 1e-9;
        end % function
        
        
        function [retval] = line_segments_intersect(...
                ax, ay, bx, by, cx, cy, dx, dy)
            % Return true if line segment ab intersections cd. NOTE: does
            % not handle degenerate cases of colinearity or colocation.
            %
            % see:
            %   http://compgeom.cs.uiuc.edu/~jeffe/teaching/373/notes/x06-sweepline.pdf
            retval = (...
                ortrace.ccw(ax, ay, cx, cy, dx, dy) ~= ...
                ortrace.ccw(bx, by, cx, cy, dx, dy) &&...
                ortrace.ccw(ax, ay, bx, by, cx, cy) ~= ...
                ortrace.ccw(ax, ay, bx, by, dx, dy));
        end % function
        
        
        function [retval] = line_segment_same_halfspace(...
                ax, ay, bx, by, cx, cy, dx, dy)
            % Returns true if line segment cd lies on the same side of the
            % line through ab. NOTE: does not handle degenerate cases of
            % colinearity or colocation.
            %
            % see:
            %   http://compgeom.cs.uiuc.edu/~jeffe/teaching/373/notes/x06-sweepline.pdf
            retval = ...
                ortrace.ccw(ax, ay, bx, by, cx, cy) == ...
                ortrace.ccw(ax, ay, bx, by, dx, dy);
        end % function
        
        
        function [retval] = ccw(a, b, c, d, e, f)
            % Returns true if the ordered points (a, b)->(c, d)->(e, f)
            % trace out a counterclockwise path.
            %
            % see:
            %   http://compgeom.cs.uiuc.edu/~jeffe/teaching/373/notes/x05-convexhull.pdf
            retval = (f-b)*(c-a) > (d-b)*(e-a);
        end % function
        
    end % methods
    
end % classdef
